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Density-functional simulations have been performed on Nass, Nag2 and Nai42 clusters 
in order to understand the experimentally observed melting properties [M. Schmidt et ai, 
Nature (London) 393, 238 (1998)]. The calculated melting temperatures are in excellent 
agreement with the experimental ones. The calculations reveal a rather subtle interplay 
between geometric and electronic shell effects, and bring out the fact that the quantum 
mechanical description of the metallic bonding is crucial for understanding quantitatively 
the variation in melting temperatures observed experimentally. 
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I. INTRODUCTION 



There is currently great interest in the statistical mechanics of finite-sized versions of systems 
that display phase transitions in the infinite-size limit. Recently, atomic clusters have provided an 
example of a finite system whose caloric curve may be measured, yielding evidence for a "melting" 
transition over a broad range of temperatures. ^'^'^'"^ The clusters involved in these experiments were 
both small (containing from a few atoms to several hundreds of atoms) and free (not supported 
on a surface). Now, according to old thermodynamic models, ^'^ a small particle should melt 
at a lower temperature than the bulk because of the effect of the surface, with the reduction in 
melting temperature being roughly proportional to 1/R, where R is the radius of the particle. This 
effect has been verified quantitatively for particles of mesoscopic size (upwards of several thousand 
atoms) supported on a surface. ^'^ More recently, however, in a series of experiments on free Na 
clusters in the size range 55 to 350, Haberland and co-workers^'^'^ found that the simple 1/R 
scaling is lost; they did observe a substantial lowering (by about 30%) of the melting temperature 
compared to bulk sodium, but accompanied by rather large size-dependent fluctuations in the 
melting temperatures. In spite of considerable theoretical effort, ^'^'^'^'^^'^^'^^ the precise form of the 
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observed fluctuations remains unexplained to date. 

As is weU known, metaUic clusters (such as sodium) possess the same first few "magic" sizes 
as atomic nuclei, A?^ = 2, 8, 20, . . . corresponding to particularly stable systems in which 
the delocalized valence electrons form closed fermionic shells. Unlike nuclei, however, metallic 
clusters also contain positively charged ions, which arc much heavier than electrons and may be 
treated classically to a good approximation, leading to the possibility that geometric packing effects 
may enter into competition with electronic shell effects in determining certain properties. The 
precise pattern of melting points observed experimentally^'^ shows maxima at sizes that correspond 
neither to magic numbers of valence electrons nor to complete Mackay icosahedra of ions, thereby 
suggesting a subtle interplay between geometric and quantum electronic effects. ^'^'^'^ 

Reliable simulations to determine the melting properties of clusters are made difficult by a 
combination of two crucial requirements: the need to compute the ionic potential-energy surface 
accurately, and the need for high statistics to converge the features in the caloric curve. Good 
statistics have been obtained using a variety of parametrized interatomic potentials, ^'-'^'^ as well 
as with a treatment of the valence electrons in the extended Thomas-Fermi (ETF)^^ and related 
density-based (DB) approximations. ^^'^^ But these attempts have failed to reproduce crucial fea- 
tures of the experimental results, such as the precise pattern of observed melting temperatures. 
An interesting observation is that in all earlier simulations the melting temperature of Na55 has 
been considerably underestimated (by about 100 K lower than the experimental value). 

Clearly, what is required is a more realistic treatment of interacting electrons, in particular, one 
that incorporates correctly electronic shell effects, which arc explicitly excluded from the above- 
mentioned work using parametrized, ETF, or DB potentials. We have recently demonstrated 
the power of ab initio density-functional theory (DFT), in the Kohn-Sham (KS) formulation, 
for simulating the melting properties of small Sn^^ and Ga^^ clusters in the size range 10—20. 
Unfortunately, the computational expense of the KS approach, combined with the relatively large 
sizes of the Na data {N = 55-350) , have so far made it difficult to perform KS simulations with high 
statistics in the size range of the experiment. Recently, a KS simulation of Na55+,^'^ performed with 
limited statistics, shows encouraging results, with an estimated melting temperature between 
300 and 350 K [r^(expt) « 290 K for Na55+]. 

In this paper, we report the first ab initio KS thermodynamic calculations up to size A^ = 142, 
within the size range of the experiment. We present simulations for Najv with N = 55, 92, and 
142, each with total sampling times of the order of 2-3 ns. For each size there is a pronounced 
peak in specific-heat curve, whose position agrees with the experimental finding to better than 8%. 
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This error is also the approximate level of statistical error for the simulation times used, suggesting 
that a KS approach is capable of quantitative predictions of melting temperatures at this level of 
accuracy or better. Analysis of the structural and dynamical information in these simulations sheds 
new light on the respective roles of geometric and electronic effects in this intriguing problem. 

II. COMPUTATIONAL DETAILS 

Our simulations have been carried out using isokinetic Born-Oppcnheimcr molecular-dynamics 
with ultrasoft pscudopotcntials in the local density approximation, as implemented in the VASP 
package.^® We have carried out calculations for up to 17 temperatures for Na55 and 14 temperatures 
for Na92 and Nai42 in the temperature range 100 to 500 K. The simulation times are 90 ps for 
temperatures below T^, and between 150-350 ps for temperatures near the melting transition 
or higher. Such high statistics have to be performed in order to achieve good convergence of 
thermodynamic indicators such as the root-mean-square bond-length fluctuation (5rms-^^ In Fig- 1, 
wc show (^ims for Na55 as a function of the total time used to compute the time averages. We 
have discarded at least the first 25 ps at each temperature to allow for thermalization. This 
plot shows that, for lower temperatures (below the transition temperature, that is, T < 240 K) 
simulation times of the order of 60-70 ps are sufficient for convergence. (For a discussion of melting 
temperatures T^, see Sec. III.) For the higher temperatures (300 K and above), at least 100-300 ps 
are required to achieve proper convergence. However, near the melting transition, namely, around 
270 and 280 K, S^ms is rather poorly converged even up to about 300 ps, indicating the need for 
higher statistics at these temperatures. Further, with increase in the cluster size, these averaging 
times may increase, making the ab initio simulations of larger clusters very difficult. 

An energy cutoff of 3.6 Ry is used for the plane-wave expansion of the wavefunction, with a 
convergence in the total energy of the order of 0.0001-0.0005 eV. Increasing the energy cutoff to, 
say, about 22 Ry causes a change of only 0.31% and 0.28%, respectively, in the binding energy and 
bond-length of the Na2 dimer. We have also checked these results against those obtained by the 
projected augmented wave (PAW) method,^'^ taking only the Is^ electrons in the core and with an 
energy cutoff of about 51 Ry. The Na2 binding energy and bond-length according to the present 
DFT method differ by less than 3% from those given by the PAW method (the PAW method 
yielding smaller energies and bond-lengths). The DFT method employed for our thermodynamic 
simulations can therefore be considered to be of comparable accuracy to the (nearly) all-electron 
method. The lengths of the simulation cells for Na55, Na92, and Nai42 are 24, 28, and 31 A, 
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FIG. 1: The root-mean-square bond-length fluctuations ^rms for Nass, as a function of the total time used 
to compute the average value. 



respectively. 

We have also performed certain tests on the Nas cluster to check the convergence of the total 
energy and forces at high temperatures. In Fig. 2, we plot the histograms of potential-energy for 
the Nag cluster at 180, 220, 300, and 450 K with a minimum of 3, 4, and 5 itcrations^"*^ for obtaining 
the self-consistent field (SCF) solution of the KS equations. Clearly, all the plots overlap, indicating 
that about 4-6 SCF iterations are sufficient to obtain a good convergence of the potential-energy. 
The force convergence test has been performed on the Nag cluster at 400 K. The forces are converged 
up to about 0.0001 eV/A in just 4 SCF iterations. Hence, in all our calculations, we have used a 
minimum of 4, and in some cases up to 5 or 6 iterations for self-consistency. 

The resulting trajectory data have been used to compute the ionic specific-heat, via a multi- 
histogram fit, as well as other thermodynamic indicators. Details can be found in Ref. 15. 
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FIG. 2: Comparison of the potential-energy distributions of the Nag cluster at 180, 220, 300, and 450 K, 
computed with minimum Ngcf = 3, 4, 5, iterations'^ used for obtaining the self-consistent field solution of 
the Kohn-Sham equations. 

III. RESULTS AND DISCUSSION 

We have carried out a substantial search for the lowest-energy ionic structures. A basin-hopping 
algorithm^^ with around 10^ basin hops was employed to generate several hundred structures 
using the second-moment approximation (SMA) parametrized interatomic potential of Li et al."^^. 
About 25 of the lowest of these structures were then used as input for relaxation within the KS 
approach. In addition, several quenches were carried out for structures selected from a long high- 
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temperature ab initio MD run, typically taken at a temperature well above the melting point (say, 
for T > 400 K). In this way, we have obtained about 20-25 distinct geometries for each size. Most 
of these structures should be close to the true ground-state, although even this method may not 
yield the correct absolute ground-state structure. The resulting lowest-energy structures found for 
the three sizes are shown in Figs. 3(a)-3(c). The lowest-energy geometry of Nass is found to be a 
two-shell Mackay icosahedron (with a very slight distortion from perfect regularity), in conformity 
with previous theory^'^^ and experimental evidence.^^ Nai42 is also found to be icosahedral, with 5 
atoms missing from the outermost shell. (A three-shell Mackay icosahedron has 147 atoms.) As to 
Na92 , the structure consists of surface growth over a slightly distorted icosahedral Nass core. The 
additional 37 surface atoms are accommodated so as to maintain a roughly spherical shape. Note 
that for Na92 , the SMA basin-hopping process has an unclear convergence with respect to the exact 
configuration of these surface atoms, and so even for the SMA potential we are unlikely to have 
found the true global ground-state structure. However, since our lowest-energy SMA structure 
turns out to be quite spherical, and since Na92 has a closed-shell configuration of valence electrons, 
we believe that the overall shape of our lowest-energy KS structure is likely to be correct. 

Next we examine the specific-heat for the three clusters investigated, which are plotted in 
Figs. 4(a)-4(c). These curves feature a single dominant peak with a width varying from about 
40 K for Na55 to about 20 K for Nai42, in general agreement with the experimental observation. 
A detailed investigation of the ionic dynamics shows this peak to correspond to a "melting-like" 
transition, from a low-temperature solidlike state, in which ions vibrate about fixed points (and 
the cluster as a whole may rotate), to a higher-temperature liquidlike state with a diffusive ionic 
dynamics. As one illustration of this, we show in Figs. 4(d)-4(f) the root-mean-squarc bond- 
length fluctuation (^rms-^^ For each size there is a sharp step in 5rms that correlates closely with 
the peak in the specific-heat. The "melting" phenomenology found is qualitatively similar to that 
found in earlier studies with simpler potentials. 

The KS melting temperatures are given along with the approximate "latent heats" L in 
Table I. Following the convention in the experimental works, ^'^ we define Tm here as the maximum 
of the peak in the specific-hcat. We have also calculated a quantity 5E defined as the average 
potential-energy of the just melted cluster with respect to the ground-state structure at T = K. 
Schmidt et al?^ have inferred from the experimental caloric curve that the melting temperature is 
strongly influenced by such an energy contribution. They showed further that follows closely 
the variation in 5E as a function of the cluster size. In Table I, we have also given and L 
calculated by the SMA parametrized potential.^'^° The striking feature of the results summarized 
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FIG. 3: Ground state geometries of (a) Na55, (b) Na92, and (c) Nai42. Structure (d) is a representative 
deformed excited-state structure of Na55. 

in this table is that only the first principles KS-based calculations reproduce, qualitatively as well 
as quantitatively, the experimentally observed variation in the melting temperatures and the latent 
heat for the three clusters investigated.^^ 

A noteworthy feature of the data in Table I is that Na92 melts at a significantly lower tempera- 
ture than Nai42. This is true both for the KS data, which include electronic shell effects explicitly, 
and for the SMA, DB, and ETF^^ data, which do not. We are thus led to conclude that geometry 
plays a significant role in the melting-point depression of Na92 : Nass and Nai42 are complete, or 
close to complete, Mackay icosahedra, while Na92 has surface growth on a two-shell icosahedral 
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TABLE I: Melting temperatures, latent heat, and 6E (see text) for NaAr given by the present Kohn- 
Sham (KS) approach, the second-moment approximation (SMA) potential, and a density-based (DB) ap- 
proach. The statistical error in the KS melting temperatures is about 8%. 
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core and is less stable. However, as we shall see, the electronic structure does play a subtle role in 
the behavior of the melting temperatures of these clusters. 

As mentioned earlier, previous simulations gave a melting temperature for Nass significantly 
lower than the experimental one. Moreover, for each of the SMA,'^'^'^ ETF,^^ and DB^^ models, 
Tm(55) is found to be significantly less than r^(142) calculated within the same model, while in 
the experiment r^(55) is very slightly greater than r^(142).^^ This discrepancy is largely removed 
within the present KS model: we find that Tjn{55) ^ r^(142), within statistical error, and that both 
melting temperatures agree with experiment. This suggests quite strongly that the high melting 
point of Na55 relative to Nai42 is due to electronic shell effects, since these are the main new element 
in our KS approach not included in previous simulations. This is particularly clear in the case of 
the DB and ETF methods. In these approaches, quantum shell effects are effectively averaged 
out as a function of size by the use of a density-dependent electron kinetic-energy functional, 
and the total cluster energy follows quite closely a smooth dependence on cluster size given by a 
liquid-drop-type formula.^^ In all other respects, however, such as the use of pseudopotentials and 
density-dependent exchange-correlation functional, the DB and ETF methods are similar to the 
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FIG. 4: Left panel: Normalized canonical specific heat for (a) Na55, (b) Na92, and (c) Nai42. Co = 
(3iV — 9/2)kB is the zero-temperature classical limit of the rotational plus vibrational canonical specific 
heat. Right panel: Root-mean-square bond-length fluctuation (5rms^^ for (d) Na55, (e) Na92, and (f) Nai42- 
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FIG. 5: Time-averaged coefBcient Cdef (see text) describing the degree of quadrupole deformation of 
Na55 (continuous line) and Na92 (dotted line). 

present KS method. We note that both Nass and Nai42 are close to "magic" systems (Nass and 
Nai38, respectively), but shell effects are relatively more important for smaller systems,^° which 
could lead to a relative enhancement in stability for Na55 and thus a relatively higher melting 
temperature. 

Further evidence for the role of quantum shell effects in melting may be obtained by examining 
the shape of the cluster before and after melting. In Fig. 5 we plot the deformation parameter edef 
for Na55, defined as edef = 2Qi/(Q2 + Q3), where Qi > Q2 > Q3 are the eigenvalues, in descending 
order, of the quadrupole tensor Qij = J2i Rii^ij- Here i and j run from 1 to 3, / runs over the 
number of ions, and Rn is the iih coordinate of ion / relative to the cluster center of mass. A 
spherical system (Qi = Q2 = Q3) has ejef = 1; a value edef > 1 indicates a quadrupole deformation 
of some type. From Fig. 5 we see that at low temperatures edef ~ Ij corresponding to the compact 
icosahedral ground-state, but as the cluster melts, the system acquires a quadrupole deformation 
with edef ~ 1.6. A more detailed investigation of this deformation with two-dimensional defor- 
mation parameters such as the Hill-Wheeler paramctcrs^'^ (not shown) indicates the cluster to be 
undergoing shape fluctuations around a net prolate deformation. A typical deformed liquidlike 
structure is shown in Fig. 3(d). A related phenomenon was observed earlier by Rytkonen et al.^^ 
for Na4o, except that the magic Na4o cluster underwent only an octupole deformation rather than 
a quadrupole deformation as observed here for the nonmagic Na55 . We do not observe statistically 
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significant quadrupole deformations for Na92 or Nai42, for which e^ef ~ 1 at all temperatures. 

Interestingly, simulations of Na55 carried out by us within the SMA model do not show a de- 
formation upon melting, but rather the cluster remains essentially spherical at all temperatures, 
^def ~ 1- Since the SMA model explicitly excludes quantum shell effects, we believe that the de- 
formation of Na55 is due to the quanta! Jahn-Teller distortion of the open-shell system of valence 
electrons. As further support for this explanation, we note that the mean prolate deformation in 
the liquidlike state agrees well with that found in jellium-model calculations for neutral Na55 that 
we have carried out, which yield a uniaxial prolate deformation with a major- to minor-axis ratio 
of i?maj/-Rmin ~ ^/cdef ~ 1-3, in close agreement with the liquidlike simulations. Evidently, the 
compact ground-state structure is favored by the possibility of geometric packing into an icosa- 
hedron, while in the nonrigid liquidlike state the cluster can lower its free energy by undergoing 
a spontaneous shape deformation. On this reasoning, the magic Na92 cluster would not be ex- 
pected to deform upon melting, consistent with the observation. On the other hand, the Nai42 
cluster would be expected to deform, at least slightly; presumably, the Jahn-Teller forces here are 
sufficiently weak that any deformation is not statistically significant. 

Note that, as mentioned earlier, the melting temperature is strongly influenced by the potential- 
energy difference SE between liquidlike and solidlike states. Therefore, even if the ground-state is 
quite spherical, as for Na55 , the melting temperature may still be influenced by important quantal 
deformation effects entering only in the liquidlike state. The KS simulations undertaken here 
incorporate all these various effects correctly. 

IV. CONCLUSION 

In conclusion, the KS approach appears to be capable of making quantitative predictions of 
melting temperatures and latent heats in Na clusters. Na55, Na92, and Nai42 are each magic or 
nearly magic, but only Na55 and Nai42 are also close to icosahedral shell closures. The fact that 
Na92 melts at a significantly lower temperature than the other two shows that geometric effects 
are very important in determining the pattern of melting temperatures observed experimentally. 
However, electronic shell effects can play an important role too, both in influencing overall binding 
energies and bond-lengths as a function of size, and indirectly via shape deformation effects that 
may arise differently in the solidlike and liquidlike states. For an accurate treatment of the metallic 
bonding, and a quantitative prediction of melting temperatures and latent heats, it is essential to 
incorporate electronic shell effects appropriately, as is possible, for example, within the KS approach 
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used here. This is especially true for the smaller sizes of clusters. There is a size regime up to 
N 150 or so where a full KS treatment is warranted. 
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